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This paper presents both preliminary experimental and numerical studies of pressure drop and heat 
transfer characteristics of Stirling engine regenerators. A test bench is designed and manufactured for 
testing different regenerators under oscillating flow conditions, while three-dimensional (3-D) numerical 
simulations are performed to numerically characterize the pressure drop phenomena through a wound 
woven wire matrix regenerator under different porosity and flow boundary conditions. 

The test bench operating condition range is initially determined based on the performance of the 
commercial, well-known Stirling engine called WhisperGen™. This oscillating flow test bench is 
essentially a symmetrical design, which allows two regenerator samples to be tested simultaneously 
under the same inflow conditions. The oscillating flow is generated by means of a linear motor which 
moves a piston in an oscillatory motion. Both the frequency and the stroke of the piston are modified to 
achieve different test conditions. 

In the numerical study, use of a FVM (finite volume method) based CFD (computational fluid dy¬ 
namics) approach for different configurations of small volume matrices leads to a derivation of a two- 
coefficient based friction factor correlation equation, which could be later implemented in an equiva¬ 
lent porous media with a confidence for future regenerator flow and heat transfer analysis. 

© 2014 Elsevier Ltd. All rights reserved. 


1. Introduction 

Nowadays Europe faces changes, driven mainly by the economic 
crisis, globalization, demand on renewable resources, and ageing. 
Focusing on the demand on renewable resources, the priority is to 
promote a sustainable growth based on efficient technologies that 
can operate with different green energy sources; one of the most 
promising technologies is the Stirling engine [1,2]. The Stirling 
engines have suitable applications including conversion of solar 
energy, co-generation, submarine and space applications. 
Regarding micro-CHP (combined heat and power), the Stirling 
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engine WhisperGen™ is an example of high efficiency with very 
low levels of chemical and noise pollution. 

The regenerator is usually the Stirling engine component that 
handles more power [3 . Therefore, the function of the regenerator 
is the key to the efficient performance of Stirling engines. Most of 
the research works in this field are carried out for the optimization 
of this component [4-6 . 

The optimization process for the Stirling engine regenerators 
mainly focus on increasing the heat transfer capacity and reducing 
pressure losses. Pressure losses through the regenerator have a 
direct influence on the Stirling engine indicated power and, 
consequently, on the electrical power output. Ibrahim and Tew [3] 
pointed out that a performance survey for small (<100 We) engines 
indicates that regenerator thermal inefficiency contributes 1.5% to 
engine thermal inefficiency while pressure drop losses contribute 
about 11% engine inefficiency. Therefore, the characterization of 
these phenomena through experimental, theoretical and numerical 
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Nomenclature 

P 

pressure, Pa 



Q 

heat flux, W/m 2 

Awr 

wetted area of regenerator, A wr = 4V dr /d h , m 2 

R 

specific gas constant, J/(kg K) 

A xr 

cross-sectional area of regenerator, A xr = V dr /L r , m 2 

Sh 

heat source term, W/m 3 

CL\ 

form resistance parameter 

Tc 

compression space temperature, K 

a 2 

surface friction parameter 

T e 

expansion space temperature, K 

Q 

friction factor = Darcy's friction 

Tg 

gas temperature, K 


factor =/ D = Ap/V2p u^ x d h /L 

hi 

gas temperature at hot side, K 

C P 

specific heat of gas at constant pressure, J/(kg K) 

Tgo 

gas temperature at cold side, K 

Cr 

specific heat of regenerator material, J/(kg K) 

^ref 

reference temperature at which physical properties are 

Cv 

specific heat of gas at constant volume, J/(kg K) 


computed, K 

dh 

matrix hydraulic diameter, m 

Twr 

regenerator matrix temperature, K 

d w 

wire diameter, m 

Twrl 

regenerator matrix temperature at hot side, K 

E 

energy per unit of mass, J/kg 

T W rO 

regenerator matrix temperature at cold side, K 

h 

convective heat transfer coefficient, W/(m 2 K) 

t 

time, s 

I< 

permeability, m 2 

U 

gas velocity, m/s 

k 

thermal conductivity of the gas, W/(m K) 

u 

mean gas velocity = rh/(pA xr ), m/s 

/<r 

thermal conductivity of the regenerator material, W/ 

14 mx 

interstitial gas velocity inside the matrix, m/s 


(m K) 

v dl - 

regenerator dead volume, V dr = V Y %, m 3 

U 

regenerator length, L r = V dr /A xr , m 

v r 

regenerator total volume, m 3 

m 

mass flow rate, kg/s 

X 

longitudinal coordinate along the regenerator length, 

Fo 

Fourier number = 4 a T /(u d h ) = 4a r /(d hv /RT ref ) 


m 

Ma 

Mach number = u/^RT ref = m^RT re f/(pA xr ) 

a Y 

thermal diffusivity of regenerator material, m 2 /s 

Pr 

Prandtl number = /i c p /k 

e 

regenerator thermal efficiency 

Re 

Reynolds number = p ud h /p = m d h //i A xr 

y 

adiabatic coefficient 

St 

Stanton number = h/(p uc p ) = M xr /(rhc p ) 

p 

viscosity of gas, Pa s 

TCR 

thermal capacity 

1v 

regenerator matrix volumetric porosity 


ratio = p r c r /(pc p ) = (p r c r T ref /p) • (y - l)/r 

p 

density of gas, kg/m 3 

n s 

engine frequency, rev/s 

Pr 

density of regenerator material, kg/m 3 


studies is crucial in co-generation applications to maximize the 
ratio between generated electrical power and source energy input. 

Most of the research works mainly focus on the empirical 
characterization of the pressure losses for stacked woven wire 
regenerator matrices. Frequently, the phenomenon is modelled as a 
case of internal flow through a conduct, and the pressure drop is 
computed by means of a friction factor correlation based on the 
Hagen-Poiseuille's law. The correlation by Kays and London [7] is 
probably the most widely used. It assumes steady, incompressible 
flow and provides a friction factor, C f correlation obtained from 
experiments performed in a perfectly pressed stacked woven wire 
screen matrices. Urieli and Berchowitz [8 proposed an equation 
based on Kay's and London's data. Seume and Simon [9] provided a 
review of the friction factor correlations for steady-state flow, and 
they also studied the compressibility effects and the characteristics 
of the oscillating flow in Stirling engine regenerators. Sodre and 
Parise [10] carried out experiments to determine the pressure drop 
through an annular conduit filled with a plain square wire-mesh 
woven-screen matrix. A corrected Ergun equation was used here 
to correlate the experimental data, considering the wall effects. 

Several experimental test benches are used to characterize both 
pressure losses and convective heat transfer in regenerators under 
unidirectional, steady flow conditions 11 . However, testing facil¬ 
ities designed for oscillatory flow are more suitable for the opera¬ 
tion conditions of Stirling regenerators. Wood et al. [12] used two 
different test benches, one for oscillating flow and other for steady 
unidirectional flow. Both test benches are designed for conven¬ 
tional heat exchangers, but they could be adapted for regenerators. 
The test bench used by Ibrahim et al. 13] consisted of an engine 
with two horizontal opposite pistons, crankshafts and individual 
drives which provide oscillating flow. Miyabe et al. [14], Tanaka 
et al. [6 , Gedeon and Wood [15] also used experimental devices 


under oscillating inflow conditions. The latter three studies [6,15] 
are classical references concerning friction factor and heat trans¬ 
fer experimental researches for stacked woven wire screens with 
variety of materials and geometries. 

Experimental works are complemented with analytical and 
numerical approaches, usually related to the flow through wire 
screen matrices. The FVM (finite volume method) appears to be 
promising numerical discretization technique for the solution of 
governing partial differential equations as indicated by Riihlich and 
Quack [5 , Gedeon and Wood [15], Ibrahim et al. [16], Tew 17] and 
others. These numerical studies suggest that the flow simulation is 
worthwhile to understand the flow behaviour and hence to char¬ 
acterize fluid flow friction for systems of potential regenerator 
applications. 

Undoubtedly, the most widely used type of regenerator is the 
metal wire matrix, mainly stainless steel. In this approach, the 
regenerator volume is filled with wires or fibers of small diameter, 
typically from 40 to 150 pm (4 x 10 _5 -15 x 10 -5 m). In the variant 
most likely used, the wires are previously woven into screen or 
meshes which are stacked adjoining to each other. There are 
different types of woven wire which lead to regenerators with 
different thermodynamic behaviour. Regenerators made of 
random-fibers or metal felts are cheaper to manufacture than 
stacked woven wire regenerators [15 . 

In the present study, a new test bench is designed for experi¬ 
mental characterization of the Stirling engine regenerator to be 
operated under oscillatory flow conditions. During the design, set¬ 
up, and preliminary tests of the test bench, a few limitations and 
measurement problems are identified. The test simplifications are 
realized to generalize the flow features and heat transfer charac¬ 
teristics however these also lead to large experimental errors. 
Therefore, a macro scale-finite volume method (FVM) based 
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numerical methodology is employed to overcome these limitations 
and/or problems to some extend and to propose more accurate 
preliminary flow solutions over a very large range of porosity and 
over a specified larger range of Reynolds number at different inflow 
velocities for a proper characterization of pressure drop phenom¬ 
ena in the Stirling engine regenerator. 

2. Main physical phenomena in Stirling engine regenerators 

The Stirling regenerator performance depends on geometrical 
parameters (d h , tv* Acr, A/vr, U V dr ), physical properties of the 
regenerator material (/< r , p r , c r ), the working fluid properties (/jl, 
c p ,y, p, k), and operating variables (rh, n s , T E , T c , p). 

The main physical phenomena to be considered in the study of 
the Stirling regenerator are summarized in the sub-sequent sec¬ 
tions with reference to research works undertaken in recent 
decades. 

2.1. Pressure drop 


occurred with the Ma in the experimental study of a non-tubular 
Stirling engine heater 18]. 

2.2. Convective heat transfer 

The heat transfer correlations are to be carefully determined as 
the friction factor correlations by taking into account the hypoth¬ 
eses under which they are obtained. Miyabe [14] carried out 
theoretical and experimental studies on heat transfer for different 
types of wire mesh. Tanaka [6] analysed the regenerator behaviour 
under oscillating flow conditions, based on experimental data for 
metallic wire meshes, sintered metals and metal foams. More 
recently, Gedeon and Wood 15] performed an extensive experi¬ 
mental study of the regenerator under oscillating flow conditions, 
analysing wire mesh and felt-type regenerators. 

A similar method to that recently used in the experimental 
study of a non-tubular Stirling engine heater [18] can also be 
applied to analyse the convective heat transfer coefficient in Stirling 
regenerators, which allows the following expression to be recom¬ 
mended as the basis for obtaining experimental correlations: 


The pressure decreases as the working gas circulates through 
the matrix of the regenerator, which in turn leads to an indicated 
power losses. The flow inside the matrix is modelled as a duct flow; 
in this case the pressure drop is obtained in analogy to 
Hagen-Poiseuille's law or to equivalent Darcy's law. Other re¬ 
searchers model the pressure losses in analogy to Forcheimer's law 
to take into account the pressure drop, not only caused by viscosity 
but also by the acceleration and deceleration effects of the fluid as it 
travels through the tortuous flow path of the porous media. 

Miyabe [14] derived generalized experimental equations of flow 
friction factor for packed wire screens with variety of geometries. 
Tanaka [6] investigated the flow characteristics of the regenerator 
in an oscillating flow for stacked woven wire and felt type matrix. 
The prediction of pressure drop loss was possible by use of the 
hydraulic diameter as the representative length defined by the 
friction factor and Reynolds number. Gedeon and Wood [15] 
derived generic correlations for friction factor based on the test 
samples of a number of stacked woven wire and felt type matrix, 
with a range of porosities. 

Following a similar method to that recently used in the exper¬ 
imental study of a non-tubular Stirling engine heater [18], theo¬ 
retically the friction factor, C f can be expressed for Stirling 
regenerators by means of the following functional relationship 
with seven dimensionless groups: 



_ n /i ^ A T* To AT 

F | Re, Ma, ~ r , H v , =-^, =X~, -= 

u h 1 wr 1 wr 1 ’ 


wr 


wr 



where T g _ and T wr _ denote, respectively, the mean temperatures of 
the gas and the regenerator matrix along the regenerator length, 
ATg = Tgi — Tgo and AT wr = T wrE — T wr q. 

The Mach number is included into Eq. (1) to account for possible 
compressibility effects, while the temperature ratios takes into 
account the spatial variations of physical properties of both the 
working gas as the regenerator material. In fact, the behaviour of 
the gas may be different during the heating and the cooling pro¬ 
cesses, but in practice it seems acceptable to obtain C f from aver¬ 
ages of measurements in both directions. 

The Eq. (1) agrees with the experimental measurements carried 
out by Su [19] for steady, compressible flow through stacked 
screens, as well as with the analyses for Stirling regenerators 
introduced by Organ [20]. However, the effects of the seven 
dimensionless group of Eq. (1) are required to be analysed with the 
experimental results because some of them could be negligible as 


wr 


St = F| Re, Pr, Ma, y, f, 1| v , 

a h 1 wr 1 wr 1 wr 


Pr c rA7g 


P ’ r h yRA7g 



2.3. Conduction heat transfer 

The convective heat transfer between the working gas and the 
regenerator matrix is related to the heat conduction inside the 
regenerator material. This relation is represented in the Eq. (2) by 
the last two non-dimensional groups, which are proportional to the 
thermal capacity ratio TCR and the Fourier number Fo, respectively. 

To study the radial heat conduction, theoretical methods are 
more common than experimental procedures [14 . 



(a) 

n 



Fig. 1. A schematic of the SRTB configuration designed. 
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Table 1 

Operating conditions of the micro-CHP Stirling engine WhisperGen™. 


Properties 

Operating conditions 


Pressure range 

Maximum pressure 
Working gas 

Heat exchanger power 

15-32 bar Hot side temperature 

40 bar Cold side temperature 

Nitrogen Maximum mass flow 

3000 W Maximum frequency 

550 °C 

80 °C 

0.02 kg/s 

25 Hz 

Table 2 

Characteristics of the SRTB heat exchangers. 

Material 

Number 
of tubes 

Inner Length 

diameter 

Maximum 

power 

Electrical heater Stainless steel 1 

Water coolers Aluminium 3 

0.01 m 0.94 m 

0.01 m 0.35 m 

5400 W 

2200 W 


For the investigated Reynolds number range of 60-350, the flow 
regime can be considered to be laminar [9] in the regenerators of 
benchmark prototypes. 


3. Experimental study 

3 A. Test bench configuration 

The SRTB (Stirling Test Bench Regenerator) is designed in 
accordance with the configuration used by Tanaka et al. [6] by 
replacing the piston drive with a linear actuator, in such a way that 
the cylinder can be completely sealed and can facilitate testing the 
flow at high pressures. The SRTB is essentially symmetrical device 
which is comprised of a central heater, two regenerators, two 
coolers and a piston as shown in Fig. 1. It is intended to obtain a 
practically steady oscillating flow in both directions of piston 
stroke. This arrangement avoids the use of high-temperature- 
resistant seals because the piston is operating at lower tempera¬ 
tures and allows the two regenerators to be tested simultaneously 
under the same flow conditions, which is important for calibration 
purposes. 

The dimensioning of SRTB is realized considering the operating 
range of the micro-CHP Stirling engine WhisperGen™ (Table 1 ). 

The linear electric actuator has a 150 mm stroke and is capable 
of exerting a continuous force of 2400 N while providing a 


maximum speed of 0.953 m/s. The main particulars of the tubular 
heat exchangers arranged in the SRTB are summarized in Table 2. 

Two regenerator specimens can be accommodated for a simul¬ 
taneous testing by means of two sample holders, formed mainly by 
a support part and a conical taper part of stainless steel (Fig. 2). The 
function of the taper part is to direct and uniform outflow of the 
heat exchangers into the matrix of the regenerator. The matrix is 
thermally isolated from the sample holder by means of mica which 
thermal conductivity is equal to 0.35 W/(m K). The gas temperature 
is measured at both input and output sections of each regenerator. 
The matrix temperature is measured at three points separated from 
each other by 120°, both in the input section and in the output of 
each regenerator. 

3.2. Data acquisition system 

In the SRTB, the following variables are monitored by means of 
transducers whose characteristics are summarized in Table 3: 7g L> 
T g o, T wrL , T wr0 , Ap, p, m and linear actuator position. 

3.3. Preliminary test results 

In the micro-CHP Stirling engine WhisperGen™, each piston has 
a stroke of 0.025 m and reaches a maximum velocity of 1.25 m/s, i.e. 
an operating frequency of 25 Hz, which is approximately 20% 
higher than the maximum frequency provided by the SRTB linear 
actuator (19 Hz). However, the frequency of the fluid oscillation is 
limited by the highest transducer response time, which corre¬ 
sponds to the MFM (mass flow meter). 

The MFM yields a direct mass flow measurement based on the 
Coriolis force principle and it is designed for measuring a one-way 
steady flow with a settling time of 0.5 s. The LA (linear actuator) 
moves in a reciprocating motion at constant speed in each direc¬ 
tion, being a triangular-wave of position vs time, so it produces a 
square-wave mass flow. During the negative half of the mass wave, 
the MFM is trying to measure a reversed mass flow, so it does not 
produce any significant measurement at all. When the square-wave 
mass flow switches from negative to positive value, the MFM starts 
to measure and it takes a time equal to the settling time of 0.5 s to 
produce a significant measurement. This settling time limits the 
maximum frequency of the LA and mass flow to 1 Hz. Because of 
this, most of the time the MFM is working in a transient regime and 
only after the settling time elapses in each positive half-wave, it 
yields valid measurements. To overcome this drawback, it is 




Fig. 2. Sample holder details. 
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Table 3 

Characteristics of the SRTB transducers. 


Transducer 

Range 

Response 

time 

Accuracy 

Quantity 

Mass flow meter 

0.4—40 g/s 

N/A 

1.0% 

1 

Pressure transducer 

0-40 bar 

5 ms 

0.2% 

2 

Differential pressure 
transducer 

-200 to 200 mbar 

5 ms 

0.1% 

2 

Thermocouple 
(Type K) 

up to 1100 °C 

50 ms 

0.5% 

16 


assumed that after the settling time elapses, there is an acquisition 
time interval in which the measurements produced by the MFM are 
inside an error band of 3% of the maximum mass flow rate (Fig. 3). 
Therefore, the valid value of the mass flow rate is assumed to be 
average of the measured values that are greater than 97% of the 
maximum values measured during the acquisition time interval. 

Due to above summarized limitations, obtaining high Re values 
to construct correlations based on Eq. (1) or (2) requires that the 
pressure be regulated adequately. The limitation is greater for Ma, 
which can only be regulated by means of the linear actuator fre¬ 
quency. Thus, in the current SRTB version, the mass flow rate 
should be calculated from indirect measurements when tests for a 
wide range of variables are necessary. 

Fig. 4 shows the variation of the instantaneous pressure drop 
test results obtained for a regenerator at gauge pressure of 16.2 bar, 
mass flow rate of 0.0086 kg/s and actuator velocity of 0.138 m/s. 
Fig. 4 also shows the instantaneous position of the linear actuator 
and the mass flow rate recorded with a settling time higher than 
the other two parameters (see Fig. 3). The time history of the 
instantaneous pressure drop changes in a triangular shape and a 
phase lag between the pressure drop and the linear actuator po¬ 
sition is observed. The measured pressure drop is found to be 
higher in the cold blow time. This result is expected because the 
pressure drop is function of the Reynolds number which becomes 
larger under constant working pressure and velocity during the 
cooling process due to increasing density. 

Although the following equations represents a very simplified 
form of the general energy conservation equation to construct a 
basis for a comprehensive study of the regenerator applications, it 
may facilitate a reasonable qualitative interpretation of the present 
experimental results here for initial design purposes. 


Settling Acquisition Direct 

Time Time Interval Mass Flow 



Fig. 5 shows a volume of the regenerator matrix of hydraulic 
diameter d h and length L r , volumetric porosity temperature of 
the solid matrix wall T wr , temperature of the working gas T g and a 
differential control volume A x dx. The solid matrix to gas heat 
transfer Q. can be expressed in terms of average heat transfer h 
coefficient for the whole regenerator volume A x L r as follows: 

Q. — fr^wr(Twr — T g ) (3) 

Starting from one-dimensional gas energy equation given by Eq. 

(4): 

tL(pE)+^[u( P E + p) + q}-S h = 0 (4) 

Considering the kinetic energy effects neglected (u 2 /2) in Eq. 
(4), and expressing p E = (c v /R)p and p E + p = c p RT g , the alterna¬ 
tive form of the energy equation is: 

jS +4 [puCpTg + q] ~ Sh = 0 (5) 

Eq. (5) is transformed into Eq. (6) assuming the mass flow is 
uniform in the matrix, the specific gas capacities are independent 
on temperature, the matrix temperature is linear and the apparent 
conductivity is not a function of x. By combining the last two as¬ 
sumptions, the term becomes as (dq/dx) = k a (d 2 T/dx 2 )~ 0. In the 
differential control volume A x dx the gas mass is the density p 
multiplied by the fraction of the control volume occupied by the 
gas, %A x dx = A xr dx and considering the solid matrix to gas heat 
transfer the gas energy equation is converted to: 

9Tg 9Tg 

PC V — + pUCp— 

The gas energy Eq. (6) is expressed in terms of regenerator 
matrix hydraulic diameter as follows: 

dT g dT g 4 , , . 

pc v —+ puc p —= 5 -/i(r wr -r g ) (7) 

For the regenerator matrix, the solid energy equation of the 
control volume is given by Eq. (8): 

dm r = hA m (T g - T wr ) ^ (8) 

The elemental regenerator matrix mass, dm r is obtained 
multiplied the regenerator material density p Y by the fraction of the 
control volume occupied by solid material. In terms the volumetric 
porosity the fraction is A xr (l - f v )dx/l v - Therefore, the regenerator 
solid energy equation is transformed into Eq. (9): 

Pr c r A xr dx (1 ~ 1v) ^ = M wr (Tg - T wr ) ^ (9) 

Hv T r 

The solid energy Eq. (9) is expressed in terms of regenerator 
matrix hydraulic diameter as follows: 

Q (l\>^ (,0) 

Fig. 6 shows a conceptual graphical representation of Eq. (6) for 
the temperature variations as function of the position and time, 
assuming that the regenerator temperature is independent of time 
at the wire surface, i.e. T wr ~T wr (x). In Fig. 6 at a time instant 
labelled as (1), the gas circulates from the hot to the cold end 
(u< 0), since dT WY /dx>0, the gas leaves heat to the regenerator 
matrix (T gx > T wr x) and the gas temperature is at a maximum value 
(0T gx /0t = 0). The gas also flows from the hot to the cold end at a 


A xr dx 


^wr(Twr — T g ) — 

Lr 


( 6 ) 
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—X— Position Pressure drop -Mass flow rate 


Fig. 4. Example of pressure drop test results. 


time labelled as (4), since T gx = T wr x and ar gx /at>0. However, it 
circulates in the opposite direction (u > 0) in the time labelled as (2) 
(Tgx — T\nyx\ ar gx /at<0) and as (3) (T gx <r wrx i dTgx/dt = 0). The 
flow direction is reversed in some intermediate time labelled as (0). 

Preliminary tests show the difficulty of measuring the temper¬ 
ature of the regenerator matrix, since it is observed that some 
thermocouples frequently give records close to the gas temperature 
values due to an insufficient contact with the matrix. In attempting 
to solve this problem it is found that the temperature of the 
regenerator matrix substantially coincides with the time average of 
the gas temperature at the corresponding end, in accordance with 
the scheme of the Fig. 6. Consequently, it seems that temperature 
measurements of the regenerator matrix are not necessary, and this 
significantly reduces the cost of the tests associated with the correct 


placement of the thermocouples. Fig. 7 shows the temperature re¬ 
sults measured in the same regenerator and conditions of Fig. 4. 

It may be noticed from Figs. 4 and 7 that the valid mass flow rate 
data are of the peak values measured during the very small acqui¬ 
sition time interval that occurs after the settling time and near to the 
reversal of the flow, which correspond approximately to the points 
labelled as (0) in Fig. 6. On the other hand, the SRTB function consists 
not only of providing facilities to obtain correlations for a given 
regenerator type, but also of performing preliminary tests to make 
comparisons among different regenerator types. For this second 
utility, the concept of regenerator thermal efficiency is of interest. 

By taking into account of the typical temperature variations of 
the oscillating flow through regenerators, it is considered adequate 
basing the thermal comparisons between regenerators on the 


Regenerator Total Volume 



Working Gas (Dead Volume) 


Regenerator Solid Matrix 



Fig. 5. Thermal analysis of regenerator. 
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Fig. 6. Conceptual scheme of flow directions and temperature distributions. 


following definitions of the heating and cooling thermal 
efficiencies: 



7gL,min ^g0,min 
■^vvrL — 7g0,min 


( 11 ) 


Table 4 

The computed efficiencies for two regenerators of 110 pm wire diameter at the 
pressure of 10 bar. 



10% 


25% 


40% 



Reg. 1 

Reg. 2 

Reg. 1 

Reg. 2 

Reg. 1 

Reg. 2 

£+ 

91.3% 

89.9% 

94.9% 

94.0% 

95.4% 

94.6% 

e~ 

99.0% 

98.1% 

99.0% 

98.7% 

99.0% 

98.3% 


Table 4 allows the comparison between the heating and cooling 
efficiencies obtained in simultaneous measurements for two re¬ 
generators of 110 pm wire diameter, at the pressure of 10 bar and 
10% of the maximum actuator velocity, which lead to a mass flow 
rate of 0.0023 kg/s at the percentile of 97%. The results obtained for 
two higher velocities (25% and 40%) are also shown. 

It is observed here that the efficiency during the heating process 
increases in both regenerators with increasing mass flow rate while 
pressure increases. However, a clear influence of the mass flow rate 
over the efficiency during the cooling process is not observed. On 
the other hand, the efficiency is higher during the cooling process 
than during the heating process in all cases, which is consistent 
with the smaller amplitude of thermal oscillations as observed in 
Fig. 7 at the regenerator cold end. This is probably due to the fact 
that the efficiency is defined in terms of temperature differences 
here. The mass flow rate is varied by controlling the gas pressure 
and velocity of the actuator. For a constant working pressure, the 
gas density is larger at the cold end of the regenerator due to lower 
temperature there. 


c - _ 7gL,max 7g0,max 

T i rp V / 

gL,max — * wrO 

where the subscripts max and min refer to the maximum and 
minimum values measured by the thermocouples at the corre¬ 
sponding sections, while T wr0 and T wrL denote the time averages of 
the regenerator matrix temperatures at each section. In order to 
diminish the noise effects on the measurements, the 97% and 3% are 
used to estimate those maximum and minimum values. 


4. Numerical study 

The lack of robustness and limitations due to flow measurement 
and testing conditions, suggest that the numerical computations 
are required for a proper characterization of pressure drop in the 
Stirling engine regenerator. The preliminary simulations are 
therefore carried out for a specified range of Reynolds number and 
volumetric porosity level in association with a simplified equivalent 
porosity media model to initially derive a two-coefficient based 
friction correlation equation to be used with a confidence for a 
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Fig. 7. Example of gas temperature measurements. 
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specified parameter range. A macro scale CFD (computational fluid 
dynamics) modelling approach is utilized with the preliminary 
objective of comparing the numerical simulation results against 
experimental test results for pressure-drop through the internal 
structure of a wound woven wire mesh matrix. By taking into 
consideration of the complexity of the flow inside the regenerator, 
the aim of the initial CFD simulations here is to study one of the 
main physical phenomena in the Stirling regenerator, the pressure 
drop. The computational modelling of the Stirling regenerator 
pressure drop losses is conducted with isothermal flow assumption 
and only one wire diameter, therefore the numbers of dimension¬ 
less group of Eq. (1) is drastically reduced. In the future, it is 
planned to extend the procedure to analyse the regenerator pres¬ 
sure drop characteristics together with the convective heat transfer. Fi S- 8 * woven wire mesh flow direction: stacked and wound. 



Flow Direction 
Stacked 


Flow Direction 
Wound 


4.1. Governing equations 


In the present numerical study, the fluid is considered to be 
incompressible and Newtonian, with constant properties, and the 
flow is assumed to be three-dimensional (3-D), laminar, isothermal 
and steady. Then the 3-D laminar form of the steady Navier-Stokes 
equations is solved for the flow field. No energy equation is solved 
for the temperature field due to isothermal flow assumption. The 
governing continuity, and momentum equations can be expressed 
in tensor notation as: 


OUj 

dxi 


0 


(13) 


dUiUj _ 1 dp d 

dXj p dX[ dXj 

where u { and p are fluid velocity and pressure, respectively while p 
and p are the density and dynamic viscosity of the fluid, 
respectively. 


dXj dXj 


-Si 


(14) 


4.2. Numerical solver 


The governing flow equation namely continuity equations and 
Navier-Strokes equations with the associated boundary conditions 
imposed are solved numerically on the basis of viscous, incom¬ 
pressible, Newtonian and laminar fluid flow assumptions in 3-D 
manner. All fluid properties (including density, viscosity, specific 
heat, and conductivity) are assumed constant. The governing 
equations are discretized and solved sequentially using a finite 
volume method (FVM) based numerical flow solver [21,22] with a 
second order upwind scheme for the discretization of the conti¬ 
nuity and momentum equations for the laminar flow solutions. The 
convergence criterion for all the velocity components and for the 
continuity is set to 10 -6 for all simulations. 


Fig. 8 illustrates the flow direction difference between a stacked 
woven wire mesh and a wound woven wire mesh configuration. 

Macroscopic measurements of the wound woven wire mesh 
matrix are conducted in order to verify the actual internal geometry 
of the matrix; it is observed that the average separation between 
mesh layers is higher than the theoretically calculated value. In 
addition, there is an expected “random” misalignment between the 
wound mesh layers as observed in Fig. 9. 

In the numerical study, the effect of porosity on the flow char¬ 
acterization in the Stirling regenerator model is investigated 
through a set of numerical simulations performed for different 
porosity level ranging from 52% to 72%, which may practically cover 
the matrix porosity range, as summarized in Table 6, with the same 
wire diameter and mesh geometry. 

Since modelling the complete regenerator matrix is too costly in 
terms of equipment and computational time, a small matrix 
domain is generated as a representative of the whole Stirling 
regenerator engine in the early stage of the numerical procedure. 
The generated 3-D geometry domain is composed of five layers 
overlapping wound mesh as shown in Fig. 10, and the entrance area 
is function of the volumetric porosity of each simulation case and 
the length of mesh is 0.001 m. 

As explained before, the aim is to study the behaviour of the 
working gas flow through the control volume. Therefore, the con¬ 
trol volume of simulation case is generated to include the fluid 
volume surrounding the wires (not the wire material) with an 
associated input distance of 0.0005 m and an output distance of 
0.0015 m to avoid backflow and to ensure that the flow is fully 
developed (Fig. 11(a) and (b)). The control volume mesh is 
composed of a hybrid mesh system, in which hexahedron mesh 
cells are utilized in the close vicinity of the entrance and exit re¬ 
gions, while tetrahedron mesh cells are utilized around the matrix 
of the mesh, with approximately 1,200,000 cells. For each simula¬ 
tion case, the maximum aspect ratio is below 2.9 and the maximum 
skewness is of 0.86. 


4.3. Model geometry and computational domain 

The numerical models are developed taking as a reference the 
Stirling regenerator matrix of the WhisperGen™ micro-CFlP Stirling 
engine. Table 5 summarizes the main features of wound woven 
wire mesh regenerator matrix used. 


Table 5 

Main mesh specifications. 


Wire 

diameter 

Mesh 

opening 

Mesh 

thickness 

Misalignment 

Porosity 

80 pm 

140 pm 

160 pm 

Random 

-60% 


4.4. Boundary conditions 

All model simulations are carried out by considering the 
following below boundary conditions as illustrated in Fig. 12: 

1. Inflow boundary: Velocity inlet boundary condition is used to 
define the uniform fluid velocity profile at the inflow boundary. 

2. Outlet boundary: Pressure outlet boundary conditions are used 
to define the static (gauge) pressure at the outlet boundary and 
to eliminate reverse flow problem. 

3. Side boundaries: Symmetry boundary conditions are used. Free- 
slip boundary conditions at the four side boundaries of the 
computational domain are imposed. The normal velocity 
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Fig. 9. Macroscopic picture of the wound woven wire mesh matrix. 


components and the normal gradients of all velocity compo¬ 
nents are assumed to have a zero value. 

4. Wall surfaces: No-slip wall boundary conditions together with 
the standard wall functions are used to define the interior wall 
boundaries between wires and flow for turbulent simulation 
cases. 

In the present preliminary computations, the CFD methodology 
adopted is restricted to isothermal, steady flow performance, 
covering the Reynolds number range representative of the instan¬ 
taneous values expected during the cyclic engine operation and one 
wound wire mesh (one diameter). 

4.5. Equivalent porous media model 

Once the initial simulations are performed for the whole ge¬ 
ometry of the wire mesh matrix, i.e. full scale computations, an 
equivalent porous media model is incorporated in the present 
numerical methodology to simulate the full regenerator flow 
behaviour for future applications with inclusion of heat transfer 
phenomena. This full regenerator model can be based on the 
pressure drop results obtained for each full scale computation. 

The geometry assumed for the equivalent porous model differs 
from that of the local detailed models as the new local mesh 
domain is modelled as an equivalent porous material. The rest of 
the control volume remains the same and the boundary conditions 
and simulation options are not modified. 

Porous media can be modelled in the present fluid flow solver 
[22] with inclusion of a momentum source term to the standard 
fluid flow equations. The source term can be composed of two parts 
namely a viscous loss term, as in Darcy's law, and an inertial loss 
term, as in Forcheimer's law. This second term dominates at high 
velocities. Laminar, transient and turbulent flow are observed in 
regenerators of benchmark Stirling engines [9 . However, at first 
approach, the flow through the regenerator wire mesh matrix is 
expected to be laminar, and Darcy's law is used to determine the 
equivalent porous media model, according to the following 
equation: 

Ap =-^u Ax (15) 

The values of permeability coefficient I< is required to be define 
the equivalent porous media model are obtained from pressure 
drop values obtained from the numerical simulations performed at 


Table 6 

Porosity values assumed in numerical models. 


Model number 

1 

2 

3 

4 

Porosity 

52% 

60% 

64% 

72% 


different velocities for the local detailed mesh models, wherein the 
thickness of the porous medium, Ax has a value of 1 mm. 

The numerical results can be used to represent a dimensionless 
chart of friction factor as a function of Reynolds number which is 
defined considering in terms of the matrix flow velocity and the 
matrix hydraulic diameter as follows: 


p„ P ^mx^h 

tvCmx = 

(16) 

P 


U 

(17) 

A — A Hv 

dh ~ dw l-1lv 

(18) 


Moreover, Darcy-Weisbach's equation allows Darcy's friction 
factor to be introduced to express the dimensionless pressure drop 
as the following variations of the Euler number referred to the 
matrix flow velocity: 

1o AP 2 =C f T (19) 

V2 p u^ x d h 

In short, it is expected that the equivalent porous media model 
lead to numerical results in agreement with the following expres¬ 
sion, which is a particular case of the Eq. (1) for unidirectional, 
incompressible, and isothermal steady flow conditions: 

C f = F(/te,1v) (20) 

The predominance of viscous or inertial effects can be assessed 
using correlations with the following form of Ergun's equation, 
where the two adjusting parameters should depend on the 
porosity: 



Fig. 10. Modelled geometry of wound woven wire mesh matrix. 
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Fig. 11. Mesh: a) wire surface mesh and b) control volume mesh. 



4.6. Numerical results 

The numerical simulations are performed for a porosity ranging 
from 52% to 72% and for Reynolds number ranging from 5 to 50 in 
the laminar flow regime to investigate the effects of porosity and 
Reynolds number on pressure drop. A two-coefficient based friction 
correlation equation is numerically obtained and its validity for the 
specified porosity and Reynolds number range is assessed in 
comparison with the experimental data available in the literature. 
The numerical results are presented in Figs. 13-17. 

Fig. 13 shows both quantitative and qualitative based pressure 
drop results obtained from the 3-D wire mesh matrix model with a 
matrix porosity level of 72%. It is clearly observed from the flow 
path lines coloured by the static pressure in Fig. 13a that the 
pressure drop occurs through the wire matrix. The pressure drop as 
a function of the length through the regenerator matrix is also 
quantitatively observed in Fig. 13b that the pressure drop behaviour 
is linear, stepwise and uniform through the regenerator matrix. In 
Fig. 13b the constant pressure drop corresponds to one obtained 
within the space between transversal woven wires and it increases 
when the flow passes through the transversal woven wires. This 
performance mainly depends on the geometrical configuration of 
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(a) 



(b) 


Fig. 13. Detailed model porosity 72% and 0.72 m/s: a) pathlines coloured by static 
pressure and b) pressure drop vs. distance. (For interpretation of the references to 
colour in this figure legend, the reader is referred to the web version of this article.) 
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the woven wire matrix. In the model of 72% volumetric porosity, the 
pressure drop is approximately 18 Pa per layer with a total pressure 
drop of 68.9 Pa (4 layers) through the model for a distance of 1 mm. 
The results obtained from the rest of the matrix configurations with 
different volumetric porosity levels also reproduce the similar 
evolution behaviour of the pressure drop. 

The plot of the pressure drop as function of the velocity through 
the different regenerator matrix configurations is also given in 
Fig. 14. It is observed here that the trend-line is approximated to a 
straight line in regards to the velocity, which corresponds to Darcy's 
law (Eq. (15)) for laminar flow in the velocity range investigated. It 
is observed an inverse influence of the volumetric porosity over the 
pressure drop, this means that the highest pressure drop is ob¬ 
tained for the volumetric porosity of 52% while the lowest pressure 
drop is obtained for the volumetric porosity of 72%. These results 
are expected considering that the volumetric porosity reduces the 
maximum velocity through the porous media (Eq. (17)) while 
increasing the hydraulic diameter (Eq. (18)), which in turn di¬ 
minishes the pressure drop (Eq. (19)). However, the influence of the 
volumetric porosity over the friction factor needs to be further 
analysed. 

These numerical results can be extrapolated to obtain the 
permeability coefficients for the equivalent porous media defined 
by Eq. (15). It should be noted here that all numerical results are 
obtained under constant working fluid properties of density and 
viscosity. 

Fig. 15 shows the path lines obtained from the 72% volumetric 
porosity equivalent porous media model at an inlet velocity of 
0.72 m/s, where the colour scale indicates the pressure drop 
through the matrix. The results for the equivalent porous media 
model based on the numerically obtained permeability coefficients 
show a good agreement with the results of the 3-D detailed wire 
mesh matrix model (see Fig. 13a), with a deviation of less than 2% 
for each model studied. 

The combined effects of volumetric porosity and the Reynolds 
number on pressure drop are illustrated in Fig. 16 in which the 
pressure drop variation is evaluated in terms of the inflow velocity 
and volumetric porosity. This figure is a 3-D representation of the 
numerical data presented in Fig. 14, where the influence of the 
volumetric porosity and the velocity over the pressure drop is 
observed. 
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Fig. 15. Equivalent model porosity 72% and 0.72 m/s: pathlines colored by static 
pressure. (For interpretation of the references to colour in this figure legend, the reader 
is referred to the web version of this article.) 


Fig. 17 shows the friction factors, which are represented as a 
function of Reynolds number for four volumetric porosities studied. 
The Reynolds number is calculated based on the hydraulic diameter 
(Eq. (16 ) while the friction factor is retrieved from the obtained 
pressure drop data (Eq. (19)). In the range of Reynolds number 
investigated, it is observed that the friction factor obtained from the 
volumetric porosity of 72% is the highest, followed by the matrix 
with a volumetric porosity of 64%. On the other hand, there are no 
significant variations in the friction factor values obtained for 52% 
and 60% volumetric porosities. Therefore, based on the present 
results it cannot be affirmed that there is a direct influence of the 
volumetric porosity over the friction factor. However, as previously 
shown in Fig. 14, there is a clear influence of the volumetric porosity 
on the pressure drop. Therefore, this initial parametric study should 
be extended to a wide range of volumetric porosities and different 


▲ Porosity 52% ♦ Porosity 60% • Porosity 64% ■ Porosity 72% 
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Fig. 14. Pressure drop versus matrix velocity at different volumetric porosity. 
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Fig. 16. Pressure drop versus matrix velocity and volumetric porosity. 


woven wire mesh configurations in order to further analyse the 
influence of volumetric porosity over the pressure drop and friction 
factor correlation equation. 

Finally, a preliminary two-coefficient based friction factor cor¬ 
relation equation is numerically derived in accordance with the 
Ergun equation (Eq. (21)) in terms of the Reynolds number. This 
equation, which is also utilized in Fig. 17 in accordance with the 
functional relation given in Eq. (21), can be defined as follows: 

230 

C f = — + 2.6 (22) 

Though the present numerical methodology in association with 
the equivalent porous media model does not completely represent 
the regenerator characterization for a very large range of para¬ 
metric values for regenerator characterization, it gives a good 
indication of determining pressure drop in the regenerator and 
hence providing a friction correlation equation to be used in con¬ 
fidence for a specified range of Reynolds number and volumetric 
porosity level. 

5. Conclusions 

In the present study, a Stirling Regenerator Test Bench (SRTB) is 
designed and constructed for Stirling engine regenerator charac¬ 
terization and a preliminary numerical procedure is presented to 
characterize pressure drop through the wound woven wire mesh 


regenerator matrices. The SRTB is designed for oscillating flow 
analysis in accordance with the configuration used by Tanaka et al. 
[6] and the operating range of the micro-CHP Stirling engine 
WhisperGen™. 

At the initial stage of the present study, the relevant geometrical 
parameters, physical properties of the regenerator material and of 
the working gas, and operating variables are identified. With regard 
to pressure drop and convective heat transfer testing, the influ¬ 
encing parameters and variables are included in a functional rela¬ 
tionship to be used as basis for obtaining experimental friction 
factor and Stanton number correlations. Some drawbacks are 
encountered during the design and setting-up of the SRTB, i.e. the 
oscillating frequency limitation, the mass flow rate measurements, 
the regenerator matrix temperature, etc. Therefore, some modifi¬ 
cation and methods are proposed during the project. 

The temperature of the regenerator matrix substantially co¬ 
incides with the time average of the gas temperature at the cor¬ 
responding end. Consequently, it seems that temperature 
measurements of the regenerator matrix are not necessary, which 
significantly reduces the cost of the tests associated with the cor¬ 
rect placement of the thermocouples. 

The concept of regenerator thermal efficiency is interesting to 
carry out comparative tests with different regenerators. The first 
tests evidence that the efficiency increases with increasing mass 
flow rate, and that the efficiency is higher for the cooling process 
than for the heating process in all cases. 
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Fig. 17. Darcy friction factor versus Reynolds number. 
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In the numerical computations on the other hand detailed wire 
mesh matrix models, i.e. the whole Stirling mesh matrix models are 
generated with different level of volumetric porosity using an 
equivalent porous media model. The initial numerical results un¬ 
derline the importance of obtaining well-characterized experi¬ 
mental data for further studies concerning heat transfer 
phenomena and the use of computational methods in the name of 
deriving a universal friction factor correlation equation for a proper 
characterization of Stirling regenerators. Though the proposed 
numerical methodology reproduces the results which are in good 
correspondence with the results of Darcy's law, the model validity 
is required to be further assessed over a very large parametric range 
and woven wire mesh configurations for future Stirling regenerator 
applications. The present numerical flow methodology also pro¬ 
vides a cost effective alternative to characterize the pressure drop 
through a woven wire mesh Stirling regenerator obtaining the 
permeability and inertial coefficient requested in a porous media 
model for further applications. 
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